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We study numerically the strong coupling limit of lattice QCD with one flavor of massless stag- 
gered quarks. We determine the complete phase diagram as a function of temperature and chemical 
potential, including a tricritical point. We clarify the nature of the low temperature dense phase, 
which is strongly bound "nuclear" matter. This strong binding is explained by the nuclear potential, 
which we measure. Finally, we determine, from this first-principle limiting case of QCD, the masses 
of "atomic nuclei" up to A = 12 "carbon" . 
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It has been a long-standing goal to determine the prop- 
erties of nuclear matter from first principles, using the 
theory which describes the interactions inside each nu- 
cleon, namely quantum chromo dynamics (QCD). Lattice 
QCD simulations have been extremely successful at de- 
termining the properties of a single proton, neutron or 
other hadron [lj, and should in principle be appropri- 
ate for studying nuclear matter as well. Great progress 
has been made in this area, with the study of nuclear 
scattering lengths and potentials 0-01, and of the ener- 
gies of two- and three-baryon systems Nevertheless, 
in these pioneering studies the quark masses are still far 
from their real-world values and the matter density is 
still very small. The ab initio study of real-world nu- 
clear matter remains a distant goal. The technical issue 
of simulating quarks which are as light as in nature is 
being resolved but two more obstacles stand in the 
way. (i) There is a large scale separation O(50) between 
the nucleons' masses and their binding energy. Accuracy 
on the latter requires excellent control over errors, (ii) 
There is a severe "sign problem" : the functional integral 
at nonzero nuclear density has an oscillatory integrand, 
which causes large numerical cancellations and altogether 
prevents its usual interpretation as a Monte Carlo sam- 
pling probability. For these two reasons, the state-of-the- 
art approach is that of effective field theory whose 
couplings still need to be derived by matching to QCD. 

Here, we consider a limit of lattice QCD where the 
above obstacles can be overcome: the limit of infinite 
coupling. Because of asymptotic freedom, the lattice 
spacing a goes to zero with the bare coupling g as a oc 
cxp( — with [3 = 6/g 2 . In the opposite limit (3 = 

we will clearly have large lattice artifacts. Moreover, dif- 
ferent discretizations of continuum QCD, all equivalent 
at weak coupling /3 1, may behave differently. We 
choose the staggered discretization of the Dirac operator 
and study the Euclidean partition function 

Z(t7 W ) = / VUVxDx e SF , (1) 



x,v=\A 



with gauge links U in SU (3) and action 

X 

(2) 



on a four-dimensional Nf x N T lattice, with antipcri- 
odic boundary conditions in Euclidean time for the 
fermions Xi periodic otherwise. The r] x p = (— l)^p<" Xp , 
i] x i = 1 arc the usual staggered phases. The chem- 
ical potential /x and an anisotropy 7 are introduced 
by multiplying the time-like gauge links U x ± | by 
7exp(±a/i/7) in the forward and backward directions, 
respectively. The usual plaquette term which accounts 
for the gluonic action is absent here, since it is mul- 
tiplied by f3. The anisotropy 7 allows the tempera- 
ture T to be varied continuously, via T = a~ 1 7 2 /N T 
at infinite coupling @. In this Letter, we consider 
only one quark species and set its mass m q to zero. 

Our model is clearly very far from continuum QCD 
which contains two light quark flavors: our "quarks" live 
on a coarse cubic crystal, come in only one flavor and 
have no spin. This causes considerable changes in the 
symmetries of the theory, its bound states and their in- 
teractions, and its phase diagram. Our model belongs to 
the family of strongly correlated fermion systems rather 
than that of quantum field theories. Still, there is am- 
ple motivation to pursue its study. First, like continuum 
QCD, it confines colored objects, and has, for m q = 0, a 
continuous global U(l) symmetry 

X(x) -S- e ie ^ a X (x), X (x) -> x{x)e l£{x)a V x , (3) 

where e(x) = (— 1)^ x " . This is the 1-flavor version of 
chiral symmetry, and is spontaneously broken in the vac- 
uum and restored at high temperature or density. Sec- 
ond, this model has been the object of analytic mean-field 
treatment since the earliest days of lattice QCD [§] , con- 
tinuing up to now (nil - fl2l |. These approximate analytic 
predictions should be checked against numerical simula- 
tions using an exact algorithm. Finally, all the obstacles 
for large (3 mentioned earlier can be tamed by exploiting 
an alternative strategy, practical when (3 = 0. 

Instead of performing the Grassmann integral in Eq.([T]) 
and obtaining a determinant which becomes a complex 
function of the gauge links U when n ^ 0, one integrates 
analytically over the gauge links first [13| . After this step, 
only colorless objects survive and propagate from one 
lattice site to a neighboring site: mesons (pions), repre- 
sented by unoriented "dimers" joining the two sites, and 
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baryons, represented by oriented triple bonds. Because 
each site hosts N c = 3 Grassmann fields \ an d 3 fields x, 
a triality constraint arises: each site is attached to exactly 
3 dimers, or is traversed by a baryon loop. Thus, baryon 
loops are self-avoiding. After performing the Grassmann 
integration, the partition function becomes a weighted 
sum over configurations of dimers (n x ^ = 0, ..,3 dimers 
for each link x, v) and self- avoiding baryon loops C flij ]: 



{nx,i> ,C} x,i> 



where a weight w(C) = p(C)^ N ^ c > exp (3kpaN T /~/) is 
associated with each baryon loop C . Here, Ng(C) is the 
number of timelike links on the loop, k is its winding 
number in this direction and p(C) = ±1 is a geometry- 
dependent sign. Thus, the weight of a configuration can 
be negative even when p = 0, which seems much worse 
than the traditional strategy. Yet, this sign problem 
can be solved by analytically resumming configurations 
where C is a baryon loop or a self-avoiding pion loop 
made of alternating single and double dimers [14(. Af- 
ter this step, the sign problem remains mild even when 
p^O [l5j], so that lattices of size 16 3 x 4 can be simulated 
using the standard technique of rcweighting at all values 
of p, thus allowing for a reliable determination of the 
full (p, T) phase diagram. The final technical difficulty 
is to devise a Monte Carlo algorithm which preserves the 
triality constraint at each site. This is achieved by the 
worm algorithm [lj|, adapted for strong coupling SU(2) 
and U(3) lattice theories in (T^J, and readily modified 
here for SU(3). It produces global updates whose high 
efficiency does not degrade as m q — >0. 

Setting the scale. We first compared, at T = p = 
and for non-zero quark mass, the results of our ap- 
proach with those of the traditional sampling of the 
fcrmion determinant by Hybrid Monte Carlo on the 
same 8 3 x 16 lattice size [18|: there was complete agree- 
ment. Then, setting m q = 0, we extracted the baryon 
mass tub from its Euclidean correlator G(r), obtaining 
am B = -1/2 log G(t + 2o)/G(t) = 2.88(1)^ close agree- 
ment with mean-field and large-iV c |19l | predictions. 
Equating our baryon mass with the real-world proton 
mass gives a ss 0.63 fm. Using instead the A mass, which 
is perhaps more appropriate for our one-flavor model, 
gives a « 0.46 fm. Alternatively, from the p meson mass 
and the pion decay constant we obtain a 0.455 fm and 
a ft* 1.4 fm, respectively. These different values give us 
an estimate of the effective coarseness of our lattice, and 
their dispersion gives an idea of the magnitude of our sys- 
tematic errors when comparing to real-world QCD: our 
model is only a caricature of the latter. 

Phase diagram. We now turn to the phase diagram as 
a function of temperature T and quark chemical poten- 
tial p [l5j . Since we take the quark mass to vanish, the 
chiral symmetry Eq.Q is exact but spontaneously bro- 
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FIG. 1: (p, T) phase diagram of 1-flavor strong coupling QCD, 
with massless staggered fermions (iV T = 4). 



ken at small (p,T), with order parameter (ipip)- When 
p = 0, a mean- field analysis predicts symmetry restora- 
tion at aT c = 5/3, whereas the Monte Carlo study of 
}2oT ] on N T = 4 lattices, extrapolated to m q = 0, finds 
aT c = 1.41(3). Here, simulating directly at m q = 0, we 
assume critical exponents of the expected 3d 0(2) uni- 
versality class, and find aT c = 1.319(2), 1.402(3), 1.417(3), 
respectively, for N T = 2, 4, 6, indicating an A r — !>oo limit 
about 15% smaller than the mean-field prediction. At 
T = the transition is strongly first-order, as we will 
see shortly. A tricritical point separates the regimes 
of first- and second-order transitions. Using finite-size 
scaling on N T = 4 lattices, we determine its location 
to be (ap T cp,aT TC p) = (0.33(3), 0.94(7)) (see Fig. Q}. 
This should be compared with the analytic prediction 
(0.577, 0.866) of [u|. The rather large difference in ^ T cp 
underlines the 0(1/ d) accuracy of a mean- field treat- 
ment, and justifies a posteriori our Monte Carlo study. 

In spite of the resemblance of Fig. [T]to the expected de- 
confincmcnt transition in massless two-flavor QCD, here 
the two phases are both confining, with pointlike mesons 
and baryons, and so the phase transition is to dense, chi- 
rally symmetric, crystalline nuclear matter. At T = 
the baryon density jumps from to 1, a saturation value 
caused by the self-avoiding nature of the baryon loops, 
which itself originates from their fermion content. In 
physical units, this represents about 4 "nuclcons" per 
fm 3 , around 25 times the real- world value. 

An intriguing feature of this T = transition is the 
value of ^ crltlcal 5 which both mean-field and an early 
Monte Carlo study [3| find much smaller than the naive 
threshold value ms/3. However, the ergodicity of the 



simulations of [14| was questioned in [21[, which was 
found to be justified in [22|. This motivated us to re- 
determine p mt (T = 0) using an improved method in- 
spired by the "snake" algorithm 23j: When two phases 
coexist, the free energy necessary to increase by a "slice" 
L x L x a the volume occupied by dense nuclear mat- 
ter can be decomposed into L 2 elementary contributions, 
looking generically like Fig. 2(a)| where one additional 
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FIG. 2: (a) Adding a baryon to grow an additional layer of 
bulk nuclear matter. Each new baryon binds to 3 nearest 
neighbors, (b) Building a first layer of nuclear matter inside 
the hadron gas, thus creating two interfaces. Each new baryon 
binds to 2 nearest neighbors. 
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static baryon is attached to 3 neighbors. We measured 
the free energy AF/T of this elementary increment on a 
large 8 3 x 16 lattice, and obtained oAf = a/4 rit = 1-78(1), 
rather close to both mean-field predictions |9| and Monte 
Carlo extrapolations [lit, b ut much smaller than arriB- 



As already recognized in [24| , the reason that /ig <mj 



must then be the presence of a strong nuclear attraction. 

Nuclear matter. Since our baryons are point-like, there 
is no conceptual difficulty in defining the nuclear poten- 
tial Vnn (R) i unlike in the real world [111 . We measured 
Vnn using again the snake algorithm, this time extending 
little by little in Euclidean time the worldlinc of a second 
baryon at distance R from the first. The result is shown 
in Fig. |3l Aside from the hard-core repulsion, there is 
indeed a strong nearest-neighbor attraction, a slight re- 
pulsion at distance ay/2, and almost no interaction be- 
yond that distance. Thus, Vnn has similar features to 
the real-world nuclear interaction, whose properties are 
commonly ascribed to subtle competition between attrac- 
tive a exchange and repulsive uj exchange. The depth of 
the minimum ~ 120 MeV and the corresponding distance 
~0.6 fm are even quantitatively plausible. This nearest- 
neighbor attraction also explains a posteriori the value 
of ^% lt : each baryon added to the dense phase binds 
with 3 nearest neighbors, which reduces the increase in 
free energy from arris to only aims + 3Vnn(o)) ~ 1.7, 
consistent with a/x^ 1 *. 

Similarly, we can predict the T = surface tension of 
nuclear matter: in a periodic cubic box, when building 
a first slice of nuclear matter with two interfaces in the 
dilute phase, each new baryon binds with only 2 nearest- 



neighbors (Fig. 2(b) ) instead of 3 in the bulk (Fig 2(a) I, 
thus increasing its free energy by | Vnn (o) I fo r an increase 
of 2a 2 in the interface area, yielding crss 2 2 _ |VNN(a)|- 

This large interface tension ~ 200 MeV/fm 2 has an 
impact on the stability of "nuclei" of various sizes and 
shapes: for a given atomic number A, those with a shape 
close to a sphere (or a cube) will have a smaller mass. 
Using the same variant of the snake algorithm, we have 
added baryons, one by one, to form such nuclei while 
measuring the successive increments in free energy. For 
A = 2 our "deutcron" binding energy is about 120 MeV: 
the real-world binding energy of ~ 2 MeV results from 



FIG. 3: Energy of a second static baryon at distance R from 
the first, i.e., (tub + Vnn(-R)), where Vnn(-R) is the nuclear 
interaction potential. The horizontal band indicates the mass 
of an isolated baryon and corresponds to Vnn = 0. At R — 
the potential is infinitely repulsive. 



delicate cancellations which do not occur in our 1-flavor 
model, and the binding energy remains of the same mag- 
nitude as the depth of Vnn- For larger A, the resulting 



Fig. 4(a) docs indeed show increased stability for nuclei 
having square (A = 4), cubic (A = 8) or parallelepipcdic 
(A = 12) shapes. Other "isomers" with different shapes, 
studied exhaustively for A = 4 and sketched Fig. |4(b)| 
have clearly larger masses. Moreover, the average mass 
per "nucleon" is well described by the first two (bulk and 
surface tension) terms of the Weizsacker phenomenolog- 
ical formula: 



m(A)/A = /4 rit + {mn) 1/z a 2 aA- 1 / 3 , 



(5) 



where a is set equal to ^-|Vnn(«)| in the figure. The 
next higher-order terms in this formula come from isospin 
and Coulomb forces, which are both absent in our model. 
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FIG. 4: (a) Mass per nucleon of A — 1, .., 12 nuclei. For 
^4 = 3,4 all possible geometric isomers are included. The solid 
line shows the parameter- free Bethe-Weizsacker Eq.(f5]), with 

the surface tension a set to |Vnn(i)|- (b) Corresponding 
nuclear geometries in order of increasing mass. 
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FIG. 5: Energy density of the pion cloud as a function of the 
Euclidean distance to a static baryon. 



Discussion and Conclusions. An interesting aspect of 
our study is the origin of the nuclear interaction. The nu- 
clcons are point-like and self-avoiding, so that only the 
hard-core repulsion is explicit. There is no direct meson 
exchange in our crude model. In a way reminiscent of the 
Casimir effect between two neutral plates, the interaction 
proceeds by the rearrangement of the pion bath caused 
by the excluded volume of the nuclcon. This rearrange- 
ment is visible Fig. [5] for one nucleon: at a neighboring 
site, the three pion lines attached to each site have fewer 
options and orient more often along the Euclidean time, 
which increases the pion energy. In fact, the nucleon 
mass aniB ~ 2.88 can be decomposed into a bare mass 
3 — 3/4 = 2.25, which is the energy increase "inside" the 
nucleon and can be assigned to the three valence quarks, 
and an energy increase ~ 0.63 in the surrounding pion 
"cloud" . When two nucleons are next to each other, the 
latter increase is limited to 10 nearest-neighbors instead 
of 2x6, which explains the attraction between them (in 
sign and roughly in magnitude). This excluded volume 
or steric effect is thus the origin of the nuclear poten- 
tial, and ultimately of nuclear stability, in our model. 
In real QCD, the pion density is not constrained as in 
Eq.(U). Nevertheless, it is going to be high at tempera- 
tures T^m^ (2(| and one should expect the same steric 
effect to enhance nuclear attraction at such temperatures. 

To summarize, in a crude model of QCD, 1-flavor lat- 
tice staggered fcrmions at strong coupling /3 = 0, we have 
been able to obtain the (/^, T) phase diagram and derive 
nuclear interactions and nuclear masses from first princi- 
ples, uncovering a simple, but universal, steric origin of 
the nuclear interaction. This model can be improved in 
many ways. One simple modification consists of giving 
a non-zero mass to the quarks: the nuclear interaction 
will weaken as the pion mass is increased, in a way which 
can be compared with effective field theories. Less sim- 
ple but feasible improvements include introducing isospin 
with a second quark flavor, and measuring the 0(f3) cor- 



rection as done analytically in [27 
our model much closer to real QCD. 
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